#the search of extremes

rm(list=ls())

# set working directory
setwd("~/Copenhagen/Copenhagen_reforestation")

Sys.setenv(LANG = "en")

#library(magrittr)


library(dplyr)
library(ggplot2)
library(EnvStats)
library(ggsignif)
library(MatchIt)
library(cobalt)

foodgroups<-c(
  #"DIETSCORE_old_nobrmlk_HH",
  "VIT.A_old_nobrmlk_HH",
  "OTHERFRUITS_old_nobrmlk_HH",
  "DARKGREENLEAF_old_nobrmlk_HH"
  #"FLESH_old_nobrmlk_HH"
  #"EGGS_old_nobrmlk_HH",
  #"MEAT_old_nobrmlk_HH",
  #"DAIRY_old_nobrmlk_HH",
  #"LEGUMES_old_nobrmlk_HH"
)

loggy<-"bin"
matching<-"full" # c("simple","genetic")
rem10k<-"norem10k"

# tries genetic:

# try10: match on reforestation
# try11: match on reforestation and navwater
# try11_ed: match on reforestation and navwater and education
# try11_simple: match on reforestation and navwatersimplebin
# try11_simple_ed: match on reforestation and navwatersimplebin and education
# try12: match on reforestation and water
# try12_ed: match on reforestation and water and education
# try13: match on none
# try14: match on navwater
# try15: match on water

tries<-c("try12","try12_ed")

#covariates, choose if you want to include deforestation
covs_diet<-c("scale(householdmembers_HH)"
             ,"scale(totunder5_HH)"
             ,"scale(householdheadage_HH)"
             #,"scale(EDUCATION_HH)"
             ,"scale(pop2000)"
             ,"scale(slope)"
             #,"scale(elevation)"
             ,"scale(accessibility)"
             ,"scale(hansen00_30pr_5k)"
             ,"scale(nightlights)"
             #,"scale(deforperc)"
)

covs_bare<-c(#"occurrence" #replaced by navwater
             "householdmembers_HH"
             ,"totunder5_HH"
             ,"householdheadage_HH"
             #,"EDUCATION_HH"
             ,"pop2000"
             ,"slope"
             #,"elevation"
             ,"accessibility"
             #,"navwater"
             ,"hansen00_30pr_5k"
             ,"nightlights"
             #,"deforperc"
)


#code 
#20jan:trying nightlights
df1<-read.csv("~/Copenhagen/Copenhagen_reforestation/Preparation files/prepared_fagan_5k_7aug23.csv")

colnames(df1)


df1$country_year<-paste(df1$survey_country_code,df1$survey_year,sep="_")
#df1$EDUCATION_HH_removednas

df1[c("area_plantation_tot_5k","area_regrowth_tot_5k")][is.na(df1[c("area_plantation_tot_5k","area_regrowth_tot_5k")])] <- 0
df1[c("navwatersimplebin")][is.na(df1[c("navwatersimplebin")])] <- 0

unique(df1$country_year)

#group by number of plantations
df2<-df1 %>% group_by(country_year)%>%
  dplyr::summarise(
    plants=mean(area_plantation_tot_5k),
    number=sum(area_plantation_tot_5k>0)
  ) %>%
  filter(plants>0) %>%
  arrange(desc(number))


df3<-df1 %>% group_by(country_year)%>%
  dplyr::summarise(
    #navwater=mean(navwater),
    navwatersimplebin=mean(navwatersimplebin)
  )


surveys <- c(unique(df2$country_year))
#surveys <- sort(surveys)
countries<-surveys[!(surveys %in% c("SN_2012","GA_2012"))]
#there is something weird going on with Gambia
df_final<-data.frame()


#countries<-countries[10:14]
#countries<-c("RW_2014","CD_2013")
#countries<-"LB_2013"

countries<-countries[c(1:9,12)]

#countries<-c("RW_2014")

for(i in 1:length(countries)){
  #i<-1
  print(i)
  country<-countries[i]
  print(country)
  # if(country %in% c("BJ_2011","CI_2011","GH_2008","GN_2012","NG_2013","TG_2013")){
  #   allresponsez<-c("LIVINGSTANDARDS_HH_removednas"
  #                   ,"MPI_HH_removednas"
  #                   ,foodgroups
  #                   )
  # } else{
  #   allresponsez<-c("LIVINGSTANDARDS_HH_removednas",
  #                   "MPI_HH_removednas")
  #   
  # }
  for(tryz in tries){
    #tryz<-"try12_ed"
    print(tryz)
    if(tryz %in% c("try10","try11","try11_simple","try11_ed","try11_simple_ed","try12", "try12_ed")){
      covs_diet2<-c(covs_diet,"scale(area_regrowth_tot_5k)")
      covs_bare2<-c(covs_bare,"area_regrowth_tot_5k")
      
      
      if(tryz == "try10"){
        covs_diet3<-covs_diet2
        covs_bare3<-covs_bare2
      } else if(tryz=="try11"){
        covs_diet3<-c(covs_diet2,"scale(navwater)")
        covs_bare3<-c(covs_bare2,"navwater")
      } else if(tryz=="try11_simple"){
        covs_diet3<-c(covs_diet2,"as.factor(navwatersimplebin)")
        covs_bare3<-c(covs_bare2,"navwatersimplebin")
      } else if(tryz=="try11_ed"){
        covs_diet3<-c(covs_diet2,"scale(navwater)","scale(EDUCATION_HH)")
        covs_bare3<-c(covs_bare2,"navwater","EDUCATION_HH")
      } else if(tryz=="try11_simple_ed"){
        covs_diet3<-c(covs_diet2,"as.factor(navwatersimplebin)","scale(EDUCATION_HH)")
        covs_bare3<-c(covs_bare2,"navwatersimplebin","EDUCATION_HH")
      } else if(tryz %in% c("try12")){
        covs_diet3<-c(covs_diet2,"scale(occurrence)")
        covs_bare3<-c(covs_bare2,"occurrence")
      } else if(tryz %in% c("try12_ed")){
        covs_diet3<-c(covs_diet2,"scale(occurrence)","scale(EDUCATION_HH)")
        covs_bare3<-c(covs_bare2,"occurrence","EDUCATION_HH")
      } else{
        print("something wrong with covss")
      }
    } else if(tryz %in% c("try13","try14","try15")){
      if(tryz == "try13"){
        covs_diet3<-covs_diet
        covs_bare3<-covs_bare
      } else if(tryz=="try14"){
        covs_diet3<-c(covs_diet,"scale(navwater)")
        covs_bare3<-c(covs_bare,"navwater")
      } else if(tryz=="try15"){
        covs_diet3<-c(covs_diet,"scale(occurrence)")
        covs_bare3<-c(covs_bare,"occurrence")
      } else{
        print("something wrong with covss")
      }
    } else{
      print("something wrong with covss")
    } 
    print(paste0("covs: ",covs_diet3))
    
    allresponsez<-c("LIVINGSTANDARDS_HH"
                    #,"MPI_HH"
                    #foodgroups
    )
    
    for (responsez in allresponsez){
      #responsez<-"LIVINGSTANDARDS_HH"
      #print(countries[i])
      print(responsez)
      #read file
      df1<-read.csv("~/Copenhagen/Copenhagen_reforestation/Preparation files/prepared_fagan_5k_7aug23.csv")
      df1$country_year<-paste(df1$survey_country_code,df1$survey_year,sep="_")
      
      df1[c("area_plantation_tot_5k","area_regrowth_tot_5k")][is.na(df1[c("area_plantation_tot_5k","area_regrowth_tot_5k")])] <- 0
      
      
      if(rem10k=="rem10k"){
        print("remove 5-10k area")
        df_10k<-read.csv("~/Copenhagen/Copenhagen_reforestation/DHS fagan overlaps/dhs10k_fagan_3aug23.csv")
        
        names(df_10k)[names(df_10k) == "area_plantation_tot"] <- "area_plantation_tot_10k"
        names(df_10k)[names(df_10k) == "area_regrowth_tot"] <- "area_regrowth_tot_10k"
        df_10k$DHSYEAR <- ifelse(df_10k$DHSCC %in% c("CI","BJ") & df_10k$DHSYEAR=="2012",2011,df_10k$DHSYEAR)
        df_10k$DHSCLUST<-as.integer(df_10k$DHSCLUST)
        df_10k$DHSYEAR<-as.integer(df_10k$DHSYEAR)
        
        df_10k_r<-df_10k %>% dplyr::select(DHSCLUST,DHSCC,DHSYEAR,area_regrowth_tot_10k,area_plantation_tot_10k)
        str(df_10k_r)
        
        class(df1$cluster_pr)
        class(df1$survey_country_code)
        class(df1$cluster_pr)
        
        df1<-merge(df1,df_10k_r,by.x=c("cluster_pr","survey_country_code","survey_year"),by.y = c("DHSCLUST","DHSCC","DHSYEAR"),all.x=TRUE)
        
        df1[c("area_plantation_tot_10k","area_regrowth_tot_10k")][is.na(df1[c("area_plantation_tot_10k","area_regrowth_tot_10k")])] <- 0
        
        
        df1$area_regrowth_tot_10k_bin<-ifelse(df1$area_regrowth_tot_10k>0,1,0)
        df1$area_plantation_tot_10k_bin<-ifelse(df1$area_plantation_tot_10k>0,1,0)
        
        
        df1_check<-df1 %>% dplyr::filter(!(area_plantation_tot_10k_bin==1 & area_plantation_tot_5k==0))
        df1_check2<-df1 %>% dplyr::filter(area_plantation_tot_10k_bin==1 & area_plantation_tot_5k==0)
        
        df1<-df1_check
        
      } else{
        print("don't remove 5-10k area")
      }
      
      
      
      
      
      #df1_check1<-df1
      #df1<-df1_check1
      
      #split per year, but this doesn't make sense anymore
      if(country=="all_2000"){
        df1 <- dplyr::filter(df1,df1$survey_year<2005)
      } else if(country=="all_2010"){
        df1 <- dplyr::filter(df1,df1$survey_year>2005)
      } else{
        print("not all")
        df1 <- dplyr::filter(df1,df1$country_year==country)
      }
      
      
      #remove dry areas
      df1<-dplyr::filter(df1,df1$latitude<11.5)
      df1_check2<-df1
      
      str(df1$latitude)
      
      #remove without latitude/longitude
      df1<-dplyr::filter(df1,!(df1$latitude<0.1 & df1$longitude<0.1))
      df1_check3<-df1
      
      
      df1[c("navwatersimplebin")][is.na(df1[c("navwatersimplebin")])] <- 0
      
      
      
      df1$villageid<-paste0(df1$survey_country_code,df1$survey_year,df1$cluster_pr)
      uniques_villageid<-unique(df1$villageid)
      #df1$MPI_HH_removednas<-as.numeric(sub(",", ".", df1$MPI_HH_removednas, fixed = TRUE))
      #df1$LIVINGSTANDARDS_HH_removednas<-as.numeric(sub(",", ".", df1$LIVINGSTANDARDS_HH_removednas, fixed = TRUE))
      #df1$EDUCATION_HH_removednas<-as.numeric(sub(",", ".", df1$EDUCATION_HH_removednas, fixed = TRUE))
      class(df1$LIVINGSTANDARDS_HH)
      
      #only select rural
      class(df1$urbanrural)
      df1<-dplyr::filter(df1,df1$urbanrural == "2")
      
      
      #set these to 0
      #colnames(df1)
      
      #calc deforestation percentage
      df1$deforperc<-df1$hansenchange00minus12/df1$hansen00_30pr_5k*100
      df1$deforperc[is.na(df1$deforperc)] <- 0
      #df1$deforperc<-
      
      plot(df1$hansen00_30pr_5k,df1$hansen12_30pr_5k)
      hist(df1$deforperc)
      
      if(loggy=="log"){
        print("logged")
        interest_pl<-"scale(area_plantation_tot_5k_log)"
        interest_reg<-"scale(area_regrowth_tot_5k_log)"
        
        df1$area_regrowth_tot_5k_log<-log(df1$area_regrowth_tot_5k+1)
        df1$area_plantation_tot_5k_log<-log(df1$area_plantation_tot_5k+1)
        
      }else if(loggy=="bin"){
        print("binary")
        interest_pl<-"area_plantation_tot_5k_bin"
        interest_reg<-"area_regrowth_tot_5k_bin"
        
        df1$area_regrowth_tot_5k_bin<-ifelse(df1$area_regrowth_tot_5k>0,1,0)
        df1$area_plantation_tot_5k_bin<-ifelse(df1$area_plantation_tot_5k>0,1,0)
        
      }else{
        print("not logged,not binary")
        interest_pl<-"scale(area_plantation_tot_5k)"
        interest_reg<-"scale(area_regrowth_tot_5k)"
        
      }
      
      #hist(log(df1$area_regrowth_tot_5k+1))
      #hist(df1$area_regrowth_tot_5k)
      #hist(df1$area_plantation_tot_5k)
      
      if(responsez == "DIETSCORE_old_nobrmlk_HH"){
        df1<-dplyr::filter(df1,df1$DIETSCORE_old_nobrmlk_HH>0)
        df1$season<-as.factor(df1$season)
        
        covs_diet4<-c(covs_diet3,"scale(LIVINGSTANDARDS_HH)","scale(agebowy_old_HH)","season")
        covs_bare4<-c(covs_bare3,"LIVINGSTANDARDS_HH","agebowy_old_HH","season")
        #covs_toadd<-c("scale(hansen00_30pr_5k)","scale(deforperc)","scale(LIVINGSTANDARDS_HH)")
        #covs_toadd_names<-c("for2000","deforperc","livst")
        fam<-"quasipoisson"
        
      } else if(responsez %in% c("LIVINGSTANDARDS_HH")) {
        covs_diet4<-covs_diet3
        covs_bare4<-c(covs_bare3)
        #covs_toadd<-c("scale(hansen00_30pr_5k)","scale(deforperc)")
        #covs_toadd_names<-c("for2000","deforperc")
        fam<-"gaussian"
        
        
      } else if(responsez %in% c("MPI_HH")) {
        "scale(EDUCATION_HH)"
        
        #covs_mpi <- covs[! covs %in% c("scale(EDUCATION_HH)")]
        #covs_bare_mpi <- covs_bare[! covs_bare %in% c("EDUCATION_HH")]
        
        covs_diet4<-covs_diet4
        covs_bare4<-covs_bare4
        #covs_toadd<-c("scale(hansen00_30pr_5k)","scale(deforperc)")
        #covs_toadd_names<-c("for2000","deforperc")
        fam<-"gaussian"
        
        
      } else{
        df1<-dplyr::filter(df1,df1$DIETSCORE_old_nobrmlk_HH>0)
        df1$season<-as.factor(df1$season)
        covs_diet4<-c(covs_diet4,"scale(LIVINGSTANDARDS_HH)","scale(agebowy_old_HH)","season")
        covs_bare4<-c(covs_bare4,"LIVINGSTANDARDS_HH","agebowy_old_HH","season")
        #covs_toadd<-c("scale(hansen00_30pr_5k)","scale(deforperc)","scale(LIVINGSTANDARDS_HH)")
        #covs_toadd_names<-c("for2000","deforperc","livst")
        fam<-"quasibinomial"
        
      }
      
      
      colnames(df1)
      
      
      inds<-c("LIVINGSTANDARDS_HH_removednas","SCHOOL_ATT_HH","SCHOOL_YRS_HH","MORTALITY_5Y_HH",
              "NUTRITION_HH","ELECTRICITY_HH","COOKINGFUEL_HH","TOILET_HH","WATER_HH","HOUSING_HH",
              "floor_HH","wall_HH","roof_HH","ASSETS_HH") 
      
      df22<-df1[c(covs_bare4,
                  inds,
                  "MPI_HH",
                  "wealthindex_HH",
                  "area_regrowth_tot_5k",
                  "area_plantation_tot_5k","area_plantation_tot_5k_bin",responsez,"villageid")]
      
      #NAs
      indx <- apply(df22, 2, function(x) any(is.na(x) | is.infinite(x)))
      
      
      
      goed<-df22[complete.cases(df22), ]
      fout<-df22[!complete.cases(df22), ]
      print(paste0(nrow(goed)," rows good"))
      print(paste0(nrow(fout)," rows fout"))
      
      #write.csv(goed,paste0("pl_goed_",country,"_",responsez,"_",tryz,".csv"))
      
      
      
      
      if(matching=="genetic"){
        pl_form_matchit <- reformulate(termlabels = c(covs_diet4), response = interest_pl, intercept = TRUE)
        
        pl_matchit<-matchit(pl_form_matchit,data=goed,
                            method="genetic",
                            pop.size=1000,
                            replace=TRUE,
                            discard="control",
                            verbose=TRUE)
        summary(pl_matchit,standardize=TRUE)
        
        saveRDS(matchit_r_final,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.rds"))
        
        love.plot(matchit_r_final,thresholds = c(m = .25),)
        
        png(paste0("pl_matchit_loveplot_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.png"))
        print(love.plot(matchit_r_final,stat = "mean.diffs", threshold = .25,binary = "raw",
                        stars = "raw",addl = addedvar))
        dev.off()
        
        df_matched_pl<-match.data(matchit_r_final,data=goed)
        write.csv(df_matched_pl,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_17oct23.csv"))
        
        
        
        
        
        
      } else if(matching=="subclass"){
        pl_form_matchit <- reformulate(termlabels = c(covs_diet4), response = interest_pl, intercept = TRUE)
        
        pl_matchit<-matchit(pl_form_matchit,data=goed,
                            method="subclass",
                            #replace=TRUE,
                            discard="control",
                            subclass=500,
                            verbose=TRUE)
        summary(pl_matchit,standardize=TRUE)
        
        saveRDS(matchit_r_final,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.rds"))
        
        love.plot(matchit_r_final,thresholds = c(m = .25),)
        
        png(paste0("pl_matchit_loveplot_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.png"))
        print(love.plot(matchit_r_final,stat = "mean.diffs", threshold = .25,binary = "raw",
                        stars = "raw",addl = addedvar))
        dev.off()
        
        df_matched_pl<-match.data(matchit_r_final,data=goed)
        write.csv(df_matched_pl,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_17oct23.csv"))
        
        
      } else if(matching=="full"){
        pl_form_matchit <- reformulate(termlabels = c(covs_diet4), response = interest_pl, intercept = TRUE)
        
        pl_matchit<-matchit(pl_form_matchit,data=goed,
                            method="full",
                            #replace=TRUE,
                            discard="control",
                            verbose=TRUE)
        summary(pl_matchit,standardize=TRUE)
        
        saveRDS(pl_matchit,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.rds"))
        
        love.plot(pl_matchit,thresholds = c(m = .25),)
        
        png(paste0("pl_matchit_loveplot_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.png"))
        print(love.plot(pl_matchit,stat = "mean.diffs", threshold = .25,binary = "raw",
                        stars = "raw"))
        dev.off()
        
        df_matched_pl<-match.data(pl_matchit,data=goed)
        write.csv(df_matched_pl,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_17oct23.csv"))
        
        
      } else if(matching=="nearestloop"){
        
        covariates<-covs_bare4
        
        goed$addtest <- 1
        goed$addtest2 <- 1
        
        cu.in<-goed
        
        varstotest <- c("addtest2",
                        covariates)
        
        v <- NULL
        vtot <- 100
        k <- 1
        addedvar<-"addtest"
        stringthatmadethecut <- NULL
        counter<-1
        
        
        while(k<100){
          print("going")
          for(i in 1:length(varstotest)){
            #i<-1
            print("start")
            print(paste0("varstotest:",varstotest))
            print(paste0("specific var to test:",varstotest[i]))
            addedvar <- c(addedvar,varstotest[i])
            print(paste0("added var: ",addedvar))
            covariates2 <- covariates[!(covariates %in% addedvar)]
            
            set.seed(42)
            matchit_r <- matchit(as.formula(paste("area_plantation_tot_5k_bin~",paste(c(covariates2),collapse="+")))
                                 ,data = cu.in
                                 #,distance="mahalanobis"
                                 ,method = "full"
                                 ,discard = "control"
                                 #,replace = TRUE
            )
            
            
            summary(matchit_r,standardize = TRUE,addlvariables = addedvar,data=cu.in)
            summy <- summary(matchit_r,standardize = TRUE,addlvariables = addedvar,data=cu.in)
            
            distance<-summy$sum.matched[1,3]
            
            #love.plot(matchit_r,stat = "mean.diffs", threshold = .25,binary = "raw",
            #          stars = "raw",addl = addedvar)
            
            nobal <- abs(summy$sum.matched[,3])>0.25
            
            summy$nn[,2]
            
            nobal_names <- as.data.frame(summy$sum.matched[,3][nobal])
            
            print("nobalnames")
            print(nobal_names)
            print(as.data.frame(summy$sum.matched[,3]))
            
            final <- nrow(nobal_names)
            
            finalmean <- max(abs(nobal_names[,1]))
            
            print("matched")
            print(summy$nn[4,2])
            print("unmatched")
            print(summy$nn[5,2])
            
            if(is.infinite(finalmean) | is.na(finalmean)){
              k<-102
              print("Balance reached!")
              #q$balance<-TRUE
              break
            }else{
              print("no balance yet")
            }
            
            
            #print(final)
            print(finalmean)
            v <- c(v,finalmean)
            if(i == length(varstotest)){
              value_batch <- min(v)
              print(paste0("end of batch - valuebatch: ",value_batch))
              #print(value_batch)
              lowest <- which.min(v)
              print(paste0("lowest of batch: ",varstotest[lowest]))
              v <- NULL
              if(value_batch<vtot){
                counter <- counter+1
                vtot<-value_batch
                x <- varstotest[lowest]
                print(paste0("halleluja, add to permanent addedvar: ",x))
                stringthatmadethecut<-c(stringthatmadethecut,x) 
                addedvar <- c("addtest",stringthatmadethecut)
                #cali<- c(cali,1.5)
                #setNames(c(cali_new, 1.5), c(names(cali_new), x))
                varstotest <- varstotest[-lowest]
              }else{
                k <- 101
                print(vtot)
                addedvar<-addedvar[1:counter]
                print("the end")
              }
            }else{
              print("not end of this batch yet")
              addedvar<-addedvar[1:counter]
            }
            
            
            
          }
        }
        
        varstotest
        addedvar<-addedvar[-1]
        
        
        covariates2 <- covariates[!(covariates %in% addedvar)]
        
        covariates2<-c(covariates2,"addtest")
        
        matchit_r_final <- matchit(as.formula(paste("area_plantation_tot_5k_bin~",paste(c(covariates2),collapse="+")))
                             ,data = cu.in
                             #,pop.size=1000 
                             ,method = "full"
                             ,discard = "control"
                             #,exact="north.south"
                             #,replace = TRUE
        )
        #library(WeightIt)
        
        summary(matchit_r_final,standardize = TRUE,addlvariables = addedvar,data=cu.in)
        
        summy <- summary(matchit_r_final,standardize = TRUE,addlvariables = addedvar,data=cu.in)
        
        love.plot(matchit_r_final,stat = "mean.diffs", threshold = .25,binary = "raw",
                  stars = "raw",addl = addedvar)
        
        nobal <- abs(summy$sum.matched[,3])>0.25
        
        summy$nn[,2]
        
        nobal_names <- as.data.frame(summy$sum.matched[,3][nobal])
        
        print("nobalnames")
        print(nobal_names)
        
        
        final_covars <- nrow(nobal_names)
        
        final_mean <- mean(abs(nobal_names[,1]))
        final_max <- max(abs(nobal_names[,1]))
        
        final_dist<-summy$sum.matched[1,3]
        
        print("matched")
        print(summy$nn[4,2])
        print("unmatched")
        print(summy$nn[5,2])
        
        final_matched<-summy$nn[4,2]
        final_unmatched<-summy$nn[5,2]
        
        
        
        row<-as.data.frame(rbind(
          c(final_dist,final_mean,final_max,final_covars,final_matched,final_unmatched)
        ))
        
        row
        
        saveRDS(matchit_r_final,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.rds"))
        
        love.plot(matchit_r_final,thresholds = c(m = .25),)
        
        png(paste0("pl_matchit_loveplot_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_7aug23.png"))
        print(love.plot(matchit_r_final,stat = "mean.diffs", threshold = .25,binary = "raw",
                        stars = "raw",addl = addedvar))
        dev.off()
        
        df_matched_pl<-match.data(matchit_r_final,data=goed)
        write.csv(df_matched_pl,paste0("pl_matchit_",country,"_",responsez,"_",matching,"_",tryz,"_",rem10k,"_5k_17oct23.csv"))
        
        
      } else{
        print("no matching")
      } 
      
      
      
       
    }
  }
  
}




